Skip to content

perf(base): speed up isR/ishom/ishom2 orthogonality checks - #213

Merged
taughz merged 1 commit into
rai-opensource:masterfrom
petercorke:perf/speed-up-isR-checks
Oct 13, 2026
Merged

taughz merged 1 commit into
rai-opensource:masterfrom
petercorke:perf/speed-up-isR-checks

Conversation

@petercorke

Copy link
Copy Markdown
Collaborator

Summary

  • isR (SO(2)/SO(3) orthogonality test used by SE3/SO3/SE2/SO2 construction with check=True) previously went through np.linalg.det/np.linalg.norm/np.eye, whose generic dispatch overhead dominates cost for such small (2x2/3x3) matrices. Added explicit fast paths for the 3x3 and 2x2 cases using a hand-written cofactor-expansion determinant and a squared-residual orthogonality check, falling back to the original generic implementation for any other size.
  • ishom/ishom2 bottom-row validity checks replaced all(T[-1,:] == np.array([...])) (array allocation + Python-level all() over a NumPy bool array) with direct scalar comparisons.
  • Net effect: isR is roughly 2x faster standalone, and SE3(T, check=True) construction is roughly 30-40% faster end to end.

Test plan

  • pytest tests/ — 343 passed
  • Verified pre-existing behaviour for symbolic (SymPy dtype=object) matrices with check=True is unchanged (was already unsupported before this change, for unrelated reasons - np.linalg.det/norm don't support dtype=object either)

🤖 Generated with Claude Code

np.linalg.det/norm and np.eye carry large dispatch overhead relative
to the actual work for 2x2/3x3 matrices, dominating SE3/SO3
constructor cost when check=True. Replace with explicit cofactor-
expansion determinants and a squared Frobenius residual for the 2x2
and 3x3 cases (falling back to the generic path otherwise), and
replace the bottom-row `all(... == np.array(...))` checks in
ishom/ishom2 with direct scalar comparisons.

~2x faster isR, ~30-40% faster SE3(T, check=True) construction.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@codecov-commenter

Copy link
Copy Markdown

⚠️ Please install the 'codecov app svg image' to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

❌ Patch coverage is 94.73684% with 1 line in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
spatialmath/base/transformsNd.py 94.73% 1 Missing ⚠️

📢 Thoughts on this report? Let us know!

@taughz
taughz merged commit 2829ced into rai-opensource:master Oct 13, 2026
13 checks passed
taughz pushed a commit that referenced this pull request Oct 13, 2026
## Summary
- `trnorm` (Gram-Schmidt re-orthonormalization used by e.g.
`SE3.__matmul__`/`@`) and `tr2adjoint` (adjoint matrix, used by
`SE3.Ad()`/`Twist3.Ad()`) both built small, fixed-size results
(3x3/4x4/6x6) via generic NumPy helpers -
`np.cross`/`unitvec`/`np.stack` in `trnorm`, `np.block` in `tr2adjoint`
- whose dispatch overhead (built for arbitrary shapes/broadcasting)
dominates cost at this size. Same pattern as the `isR`/`ishom` speedup
in #213.
- `trnorm`: replaced with explicit scalar cross-product/normalization
arithmetic and direct in-place array construction.
- `tr2adjoint`: replaced `np.block([[...]])` with a pre-allocated
`np.zeros((6,6), dtype=T.dtype)` and direct slice assignment.
`dtype=T.dtype` is preserved specifically so the function's documented
`:SymPy: supported` contract still holds for symbolic (`dtype=object`)
input.
- Net effect (isolated `timeit`, `trnorm`/`tr2adjoint` called directly):
`trnorm` ~12x faster (27.4us -> 2.3us), `tr2adjoint` ~2.7x faster (8.8us
-> 3.3us). End to end, `SE3 @ SE3` (which normalizes via `trnorm`) drops
from ~30us to ~4us. `SE3.Ad()`/`base.tr2adjoint()` calls on an
already-known SE3/rotation matrix land around ~3.3us.
- Note: `Twist3.Ad()` itself is not meaningfully faster from this change
- its cost is dominated by the separate `self.SE3()` -> `trexp()`
conversion step (matrix exponential), not by `tr2adjoint`. That's a
different, inherent cost (same class as `Twist3 * Twist3`'s "product of
exponentials"), not addressed here.

## Test plan
- [x] `pytest tests/` - 343 passed
- [x] Verified bit-for-bit numeric equivalence against the prior
implementation over 200 random SO(3)/SE(3) trials (both `trnorm` and
`tr2adjoint`, both 3x3 and 4x4 forms)
- [x] Verified symbolic (SymPy, `dtype=object`) equivalence for
`tr2adjoint` specifically, since it's the one with a documented `:SymPy:
supported` contract

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
taughz pushed a commit that referenced this pull request Oct 13, 2026
## Summary
- `qqmul` (Hamilton product) and `qvmul` (vector rotation) both used
`np.dot`/`np.cross` on 3-vectors, whose generic dispatch overhead -
built for arbitrary shapes/broadcasting - dominates cost at this size.
Same pattern as `isR`/`trnorm`/`tr2adjoint` in #213/#214.
- `qqmul`: replaced with explicit scalar Hamilton-product arithmetic.
Bit-identical to the prior result (verified to ~1 ULP over 2000 random
trials).
- `qvmul`: replaced the `q * pure(v) * conj(q)` sandwich (two full
Hamilton products through `qqmul`/`qpure`/`qconj`, each with its own
`getvector` re-validation, and each Hamilton product wasting work on a
zero scalar part) with the closed-form rotation identity `v' = v + 2s(w
x v) + 2 w x (w x v)` for `q = (s, w)`. This is a genuinely different
(though mathematically equivalent) formula, not just a re-expression of
the same one - so it is **not** bit-identical, verified instead to
~1e-14 absolute over 2000 random trials (with `v` scaled to magnitude
~10), well within floating-point noise.
- Net effect: `qqmul` ~10x faster (15.0us -> 1.4us isolated), `qvmul`
~24x faster (33.4us -> 1.4us isolated). End to end: `Q1 * v` drops from
~33us to ~4us (~9x). `Q1 * Q2` drops from ~19us to ~10us (~2x) - `qqmul`
itself is no longer the bottleneck there; the remainder is
`UnitQuaternion` construction overhead (in particular `qunit()`'s
`np.linalg.norm`/`np.r_` calls), a separate, un-addressed candidate for
a possible future fix in the same spirit, not attempted here.

## Test plan
- [x] `pytest tests/` - 343 passed
- [x] `qqmul` verified against the prior implementation over 2000 random
quaternion pairs (max diff ~1 ULP, i.e. bit-identical modulo
floating-point operation ordering)
- [x] `qvmul` verified against the prior implementation over 2000 random
(quaternion, vector) pairs, vector magnitude ~10 (max abs diff ~1.4e-14)

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
taughz pushed a commit that referenced this pull request Oct 13, 2026
…#216)

## Summary

`trexp`/`trlog` are much slower than `rodrigues` (~20us and ~17us vs
~7us), but unlike the previous three perf PRs (#213, #214, #215) they
don't have one dominant wasteful generic-NumPy call to fix - they're
deep call chains (`trexp` -> `isskewa` -> `vexa` -> `iszerovec` ->
`unittwist_norm` -> `rodrigues` -> `skew` -> `rt2tr` -> `ishom` ->
`isR`) where cost is spread thin across many small layers, each
independently re-validating/re-extracting. This PR pulls out the two
safe, verified wins found in that chain; it does **not** close the full
gap to `rodrigues` - see "further speedup opportunities" below.

- **`isskewa`** (the validity check `trexp` runs on its se(3) input) had
the exact same bug `ishom` had before #213: `np.linalg.norm()` on a
small fixed-size matrix, plus an `all(S[-1,:] == 0)`
array-allocation-and-compare for checking the bottom row is zero. Fixed
the same way (explicit sum-of-squares residual + scalar comparisons).
~1.9x faster standalone.
- **`np.eye(3)`** was rebuilt from scratch on every call in `rodrigues`
(1x), `trexp`'s `V`-matrix construction (1x), and `trlog` (2x) - despite
only ever being used as a read-only operand in an
addition/multiplication that produces a new array. Replaced with a
module-level constant `_EYE3` at those call sites only. **Not** used at
the four zero-motion/zero-rotation early-return `np.eye(N)` calls in
`rodrigues`/`trexp` (e.g. `if iszerovec(w): return np.eye(3)`) - those
hand the array object directly to the caller, and aliasing them to a
shared constant would risk a real bug if a caller ever mutated the
returned matrix in place. Left as fresh arrays.
- **`isskew`** (the plain, non-augmented so(n) check) was also tried
with the same explicit-arithmetic treatment as `isskewa`, but did
**not** show a reliable win under careful (min-of-repeats) benchmarking
- the only cost it avoids is `np.linalg.norm` on an already-cheap `S +
S.T`, not enough margin to reliably beat by hand-unrolling the indexing.
Reverted to the original implementation rather than keep an unproven
"fix" - same judgment call as an earlier fully-scalar `isR` variant that
was tried and discarded for the same reason.

Net effect, measured old-vs-new in the same process, min of 9 repeats
each (single-run `timeit` numbers were misleading here - see commit
message): **isskewa ~1.9x, rodrigues ~1.13x, trexp ~1.05x, trlog
~1.10x**. Real, but modest compared to #213-#215.

## Further speedup opportunities (not attempted here)

A bigger win on `trexp`/`trlog` would mean collapsing the redundant
`getvector`/re-validation layers in the hot path into a more direct,
fused implementation - e.g. `trexp(se3_matrix)` currently re-extracts
and re-validates the same data through `isskewa` -> `vexa` ->
`iszerovec` -> `unittwist_norm` independently, each with its own
`getvector` call. That's a larger, riskier change: this code has
clearly-deliberate handling for numerically delicate edge cases
(near-zero rotation, near-pi singularity, caller-supplied `theta`), so a
fused rewrite needs much more careful edge-case test coverage than the
drop-in fixes here. Not scoped or attempted in this PR.

Two smaller, unrelated candidates surfaced but out of scope for this PR:
- `np.linalg.norm` is used directly (instead of the project's own
`norm()`, which is faster for small vectors per its own docstring) in
several other places - `iszerovec`, `isunitvec`, `qnorm`/`qunit` in
`quaternions.py`, and a few spots in
`geom2d.py`/`geom3d.py`/`pose3d.py`. Worth noting: `iszerovec` was
checked directly as part of this investigation, and swapping it to the
project's `norm()` did **not** show a reliable win either (consistent
with the project's own timing script, which shows `base.norm(R6)`
measuring *slower* than `np.linalg.norm(R6)` for a 6-vector) - so this
isn't a guaranteed win either, just an inconsistency worth someone
measuring properly before touching.
- `np.eye(N)` is rebuilt from scratch in ~20 other places across the
codebase (`pose2d.py`, `pose3d.py`, `transforms2d.py`,
`transformsNd.py`, and other unrelated functions in `transforms3d.py`).
None are in the `trexp`/`trlog`/`rodrigues` path touched here; a broader
sweep would need the same operand-vs-return-value aliasing care taken in
this PR, function by function.

## Test plan
- [x] `pytest tests/` - 343 passed
- [x] `isskewa` verified against the prior implementation over 3000
randomized trials (skew/non-skew, near-boundary perturbations)
- [x] `rodrigues`, `trexp`, `trlog` (both matrix and twist-vector return
forms) verified bit-identical (max diff 0.0) against reference copies of
the prior implementations over 1000-2000 randomized trials each

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants